Comparações entre três ou mais condições dependentes e independentes
Paulo S. P Silveira (paulo.silveira@fm.usp.br)
Koichi Sameshima (koichi.sameshima@fm.usp.br)
José O. Siqueira (jose.siqueira@fm.usp.br)
v20200408.1936
Objetivos
- reconhecer e mencionar propriedades da distribuição \(F\).
- reconhecer as indicações e aplicar ANOVA para três ou mais condições independentes.
- reconhecer as indicações e aplicar ANOVA para três ou mais condições dependentes (medidas repetidas).
- definir hipóteses estatísticas nula e alternativa.
- executar e interpretar os testes estatísticos omnibus e post hoc.
- executar e interpretar estatísticas de tamanho de efeito.
Preparação
Os exemplos aqui apresentados estão disponíveis. Caso queira usá-los, crie um projeto, coloque o arquivo desta aula e os seguintes arquivos na pasta do mesmo:
- Animacao_F.R
- Testes_t_2a2.R
- HardyWeinberg.R
- Nutricao.xlsx para o teste \(t\)
- ANOVA1f_indep_FisherWhite_sodio.R
- ANOVA1f_indep_Welch_sodio.R
- Nutricao3.xlsx para ANOVA independente
- ANOVA1f_dep_balanc_sodio.R
- Nutricao3par.xlsx para ANOVA relacionada, balanceada
- ANOVA1f_dep_desbalanc_sodio.R
- Nutricao3parD.xlsx para ANOVA relacionada, desbalanceada
Sobre os métodos tradicionais
ANOVA (sigla, do inglês, Analysis of Variance) utiliza a estatística \(F\).
- A ANOVA de Fisher é uma extensão do teste \(t\), com:
- VI (variável independente) nominal com três ou mais categorias (politômica nominal, fator)
- VD (variável dependente) quantitativa (unifatorial, One way ANOVA)
- Pode ser:
- independente (fator entre participantes)
- relacionada (fator intraparticipantes)
O tamanho do efeito é expresso pelo eta ao quadrado (\(\eta^2\)).
|
Tradicionalmente exigia-se, para a ANOVA clássica, que a distribuição populacional da variável fosse normal para cada condição e que existisse homocedasticidade entre os grupos. Estas premissas eram testadas:
Quando as premissas não eram atendidas, recomendava-se o teste não paramétrico H de Kruskal-Wallis. |
ANOVA é relativamente robusta às violações das suposições de normalidade, homocedasticidade e outliers.
Se o tamanho da amostra é pequeno, ou há heterocedasticidade, desbalanceamento, ou assimetria da VD, deve-se considerar a execução de uma ANOVA de Welch, ou ANOVA não-paramétrica ou ANOVA com reamostragem (bootstrapping).
- é necessário seguir as premissas com todo o rigor?
Ilustramos com uma analogia, recordando o equilíbrio de Hardy-Weinberg.
The Hardy-Weinberg Law of Genetic Equilibrium
In 1908 G. Hardy and W. Weinberg independently proposed that the frequency of alleles and genotypes in a population will remain constant from generation to generation if the population is stable and in genetic equilibrium. Five conditions are required in order for a population to remain at Hardy-Weinberg equilibrium:
- A large breeding population
- Random mating
- No change in allelic frequency due to mutation
- No immigration or emigration
- No natural selection
Violando as condições de H-W, em simulação com:
- 500 indivíduos,
- mutação alterando a frequência dos alelos,
- selecionando para o alelo A.
source("HardyWeinberg.R")Observamos, por exemplo: 
As condições matemáticas de população muito grande (ou infinita), sem mutação e sem seleção garantem a população no equilíbrio de Hardy-Weinberg; porém, dentro de certos limites, ainda encontramos predição aceitável para o Equilíbrio de Hardy-Weinberg e poderíamos utilizá-los em populações biológicas.
- o que é independência?
Um delineamento independente é aquele em que cada unidade experimental só é submetida a uma condição experimental.
Por exemplo, três grupos com 5 indivíduos cada um (total de 15 participantes, nomeados de A a J) são alocados assim:
| Condição 1 | Condição 2 | Condição 3 | |
| A | F | K | |
| B | G | L | |
| C | H | M | |
| D | I | N | |
| E | J | O |
Em um delineamento relacionado, o mesmo indivíduo é submetido a todas as condições experimentais. Os mesmos 15 participantes seriam alocados assim:
| Condição 1 | Condição 2 | Condição 3 | |
| A | A | A | |
| B | B | B | |
| C | C | C | |
| D | D | D | |
| E | E | E | |
| F | F | F | |
| G | G | G | |
| H | H | H | |
| I | I | I | |
| J | J | J | |
| K | K | K | |
| L | L | L | |
| M | M | M | |
| N | N | N | |
| O | O | O |
- o que é balanceamento?
Em delineamentos em que os grupos são balanceados, o número de participantes em cada condição é aproximadamente igual. Por exemplo, com 15 participantes em ANOVA independente, desbalanceada, poderíamos ter algo como.
| Condição 1 | Condição 2 | Condição 3 | |
| A | H | K | |
| B | I | L | |
| C | J | M | |
| D | N | ||
| E | O | ||
| F | |||
| G |
Em um delineamento intra-participantes o desbalanceamento ocorre quando um dos participantes deixa de aparecer em uma das condições experimentais, por exemplo:
| Condição 1 | Condição 2 | Condição 3 | |
| A | A | A | |
| B | B | ||
| C | C | C | |
| D | D | D | |
| E | E | E | |
| F | F | F | |
| G | G | G | |
| H | H | H | |
| I | I | I | |
| J | J | ||
| K | K | K | |
| L | L | L | |
| M | M | M | |
| N | N | N | |
| O | O | O |
Basta uma participação faltando (missing) para caracterizar o desbalanceamento. O procedimento com as estatísticas padrão do R não funcionarão e temos que usar um modelo com efeitos aleatórios; lidaremos com esta situação adiante.
Distribuição \(F\) de Fisher-Snedecor
Esta distribuição é uma generalização da estatística \(t\) para três ou mais grupos. Descreve probabilidades que dependem da razão entre duas variâncias, e portanto considera graus de liberdade (\(\nu\), letra grega ni ) para o numerador e para o denominador.
Para a comparação de três ou mais condições, com base na estatística \(F\), os graus de liberdade do numerador dependem do número de condições e os do denominador dependem do tamanho da amostra.
Familiarize-se com a distribuição \(F\), observando Animacao_F.R:
- \(F\) é um valor maior que zero.
- sob \(H_0\) a distribuição \(F\) tem um parâmetro de não centralidade (
ncp) igual a zero; sob \(H_1\) o parâmetro de não centralidade é maior do que zero. - esta distribuição é assimétrica.
- consideramos somente a cauda superior.
- localize \(\alpha\) e \(\beta\).
- há dois valores para graus de liberdade: para o numerador (número de grupos - 1) e denominador (número de sujeitos - número de grupos).
- observe o que acontece com a distribuição sob \(H_0\) e com o valor de \(F\) crítico (\(F_c\)) à medida que os graus de liberdade aumentam.
O delineamento dos estudos, o tipo de variável e, consequentemente, a estatística adequada mudam, mas o problema é sempre o mesmo: incerteza porque lidamos com uma amostra.
A diferença, aqui, é que teremos que lidar com 3 ou mais amostras simultaneamente.
A hipótese nula é pela igualdade de todos os \(m\) grupos. caso não rejeitemos \(H_0\): 
que pelo menos uma destoe das demais (i.e., basta que uma das condições ser diferente das demais para rejeitarmos \(H_0\).), por exemplo: 
teste \(t\) (quando eram duas turmas)
Brendon Small e Coach McGuirk fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as duas turmas.
Aplicamos o teste \(t\) de Satterthwaite (a.k.a. Welch), robusto e default do R:
# dados
library(readxl)
Dtfrm <- read_excel("Nutricao.xlsx")
# significancia estatistica
t_out <- t.test(Sodium ~ Instructor, data = Dtfrm)
print(t_out)
Welch Two Sample t-test
data: Sodium by Instructor
t = 0.76722, df = 34.893, p-value = 0.4481
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-67.91132 150.41132
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk
1287.50 1246.25 testes \(t\) quando temos três ou mais turmas
Suponha que uma terceira classe junte-se ao experimento (Melissa Robins). Podemos, então, verificar as diferenças destas três condições experimentais com testes \(t\), comparando-se Brendon Small com Coach McGuirk, Brendon Small com Melissa Robins e Coach McGuirk com Melissa Robins?
Quantos testes \(t\) precisaremos de acordo com o número de grupos?
É a combinatória do número de grupos (\(m\)) dois a dois: \[{m \choose 2} = { {m!} \over {2! (m-2)!} }\]
- para 2 -grupos: \({2 \choose 2} = { {2!} \over {2! (2-2)!} } = 1~\text{teste}\)
- para 3 grupos: \({3 \choose 2} = { {3!} \over {2! (3-2)!} } = 3~\text{testes}\)
- para 4 grupos: \({4 \choose 2} = { {4!} \over {2! (4-2)!} } = 6~\text{testes}\)
- para 6 grupos: \({6 \choose 2} = { {6!} \over {2! (6-2)!} } = 15~\text{testes}\)
Graficamente vemos:
# Testes_t_2a2.R
# Ilustra o numero de testes t necessarios
# para comparar grupos par a par
combinacoes <- c()
grupos <- c()
for (g in 2:20)
{
grupos <- c(grupos,g)
combinacoes <- c(combinacoes,choose(g,2))
}
plot(grupos,combinacoes,
xlab="Número de grupos",
ylab="Número de testes t feitos par a par",
lwd=2,
type="b")
Embora o número cresça rapidamente, podemos não nos impressionar pois temos computadores para o trabalho repetitivo. Há, porém, um problema mais grave: as probabilidades dos erros do tipo I (\(\alpha\)) e do tipo II (\(\beta\)) se acumulam. Quando escolhemos \(\alpha\), a probabilidade de rejeitar incorretamente a hipótese nula, o máximo valor que \(p\) pode assumir, a probabilidade deste erro para \(m\) grupos é:
\[P(\alpha | m) = 1-(1-\alpha)^{m \choose 2}\]
Também gostaríamos de manter a probabilidade de rejeitar corretamente a hipótese nula, o poder do teste, \(1-\beta\); para \(m\) grupos é:
\[P(\beta | m) = (1-\beta) ^ {m \choose 2}\]
Graficamente podemos observar o que acontece com o número crescente de pares de teste \(t\) necessários, considerando os tradicionais \(\alpha=0.05\) e \(\beta=0.1\):
# Testes_t_AlfaPoder.R
# Ilustra o valor de alfa e 1-beta
# com numero crescente de testes
source("friendlycolor.R")
cor_alfa <- friendlycolor(8)
cor_poder <- friendlycolor(20)
alfa_escolhido <- 0.05
beta_escolhido <- 0.20
alfa <- c()
poder <- c()
grupos <- c()
for (g in 2:20)
{
grupos <- c(grupos,g)
combinacoes <- choose(g,2)
alfa <- c(alfa,1-(1-alfa_escolhido)^combinacoes)
poder <- c(poder,(1-beta_escolhido)^combinacoes)
}
plot(grupos, poder,
xlab="Número de grupos",
ylab="Probabilidade",
ylim=c(0,1),
lwd=2, lty=2, col=cor_poder,
type="b")
lines(grupos,alfa,
lwd=2, lty=1, col=cor_alfa,
type="b")
legend("right",
c("alfa", "poder"),
col=c(cor_alfa,cor_poder),
lwd=c(2,2),
lty=c(1,2),
box.lwd=0, bg="transparent")
Portanto, a probabilidade de erro do tipo I cresce rapidamente para quase \(100\%\) e o poder de seus testes combinados vai para zero. Na prática, dificilmente teremos mais do que 6 grupos mas, ainda assim, teríamos \(\alpha \approx 53.7\%\) e poder \(1-\beta \approx 3.5\%\), valores totalmente inaceitáveis para uma boa análise.
Podemos resolver com vários testes \(t\)? A resposta é: 
NÃO!
Precisamos testar tudo simultaneamente paramanter \(\alpha\) no valor pretendido e
preservar \(1-\beta\).
Análise da variância
Esta análise, que tem o acrônimo ANOVA (do inglês, Analysis of Variance) utiliza apenas variâncias, mas…
… as conclusões são sobre as médias dos \(m\) grupos.
Como?
Princípio do teste ANOVA unifatorial
Consideremos que existam três grupos com distribuição normal. De cada um deles, retiramos uma amostra: 
Não havendo problemas, as três amostras reproduzem a distribuição, média e variância das populações das quais se originaram.
A distribuição \(F\) é dada por \[F = {{S_E^2}\over{S_D^2}}\] onde \(S_E^2\) é a variância entre os grupos e \(S_D^2\) é a variância dentro dos grupos.
Variância entre (between) os grupos
A variância entre os grupos presume que as médias amostrais (\(\bar{X}_m\)) refletem as respectivas médias populacionais (\(\mu_m\)): 
Como sempre, não temos esta certeza e lidamos somente com a informação das amostras, \(\bar{X}_m\): 
Cada média é um número e, portanto, podemos calcular a variância destes números: 
Sendo assim, para \(F = {{S_E^2}\over{S_D^2}}\), a estatística \(F\) aumenta quando a variância entre os grupos (ou condições) aumentar.
Variância dentro (within) dos grupos
A variância dentro dos grupos é uma medida da variância total, desconsiderando a média de cada uma das condições. Cada amostra tem sua própria distribuição (presumivelmente, reflexo da distribuição da população de onde veio): 
Como sempre, não temos esta certeza e lidaremos somente com a informação das amostras, \(\S_m\): 
Mesclando as três distribuições, estimamos a variância dentro dos grupos, \(S_D^2\), uma medida de quanto, como um todo, a variável é dispersa: 
Caso a variância em cada condição seja maior, 
esta variância será refletida em \(S_D^2\): 
Sendo assim, para \(F = {{S_E^2}\over{S_D^2}}\), a estatística \(F\) diminui quando a variância dentro os grupos (ou condições) aumentar.
Comportamento de \(F = {{S_E^2}\over{S_D^2}}\)
É fácil imaginar o comportamento da estatística \(F\) combinando-se o que pode acontecer com \(S_E^2\) e \(S_D^2\):
Desta forma, Ronald Fisher inventou uma forma de comparar médias entre várias condições, simultaneamente, utilizando somente a comparação entre duas variâncias, com o numerador refletindo a dispersão das médias e o denominador refletindo a dispersão do fenômeno em estudo. 
Métodos Robustos
Os métodos aqui apresentados são robustos à falta de homocedasticidade. A ausência de normalidade é mais difícil de ser contornada em amostras pequenas. Podemos prescindir da suposição de normalidade para amostras maiores (tamanhos mínimos são vistos adiante) utilizando o teorema central do limite.
ANOVA para condições independentes
Utilizada quando os participantes são avaliados em somente uma das condições experimentais, i.e., um delineamento independente ou entre participantes.
- suposições
- independência entre as unidades experimentais,
- normalidade da VD em todas as condições,
- homocedasticidade da VD entre todas as condições.
|
Normalidade e balanceamento:
Então: - se o número de unidades experimentais é igual a pelo menos 12 em cada grupo, a suposição de normalidade não é necessária.
O balanceamento não precisa ser estrito:
Então:
|
|
Homocedasticidade: Em geral, quando você tiver tamanhos de amostras iguais, a suposição de heterocedasticidade não será um grande problema.
- Dancey & Reidy (2019), p. 472-3
- Johnson, R. & Wichern, D. (2007) Applied Multivariate Statistical Analysis. 6 th ed. NJ: Prentice-Hall, p. 291
- Moore, D. (1995) The basic practice of statistics. New York: W. H. Freeman and Company,
citado por Norusis, M. (2009) PASW Statistics 18: Statistics Procedures Companion. NJ: Prentice-Hall, p. 148 |
- situação
O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.
O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.
Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.
- planejamento
Brendon Small, Coach McGuirk e Melissa Robins fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as três turmas.
- hipóteses nula e alternativa
As três classes receberam diferentes programas de educação nutricional. A ingestão média de sódio é a mesma (populacionalmente) para os três programas?
\[H_0: \mu_{\text{Small}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\] \[H_1: \exists |\mu_i \ne \mu_j; i \ne j; i,j=1,2,3\]
|
Esta é uma forma matemática para escrever a hipótese alternativa, lida como “Existe pelo menos alguma média \(\mu_i\), diferente de outra média \(\mu_j\), com \(i\) e \(j\) assumindo os valores \(1\), \(2\) ou \(3\)”. É dizer que \(\mu_1 \ne \mu_2\) OU \(\mu_1 \ne \mu_3\) OU \(\mu_2 \ne \mu_3\). Aparece, por aí, como: \[H_1: \mu_1 \ne \mu_2 \ne \mu_3\] tentando indicar as diferenças, mas esta forma não funciona porque sugere que a rejeição de \(H_0\) implica em que todos os grupos sejam diferentes entre si, como \(\mu_1 \ne \mu_2\) E \(\mu_1 \ne \mu_3\) E \(\mu_2 \ne \mu_3\). Basta que uma condição tenha média estatisticamente diferente das demais para que se rejeite \(H_0\) com ANOVA. Em português: \[H_1: \text{Existe pelo menos um grupo diferente de algum outro.}\] |
Operacionalização do teste
ANOVA unifatorial independente de Fisher com ajuste para heterocedasticidade de White
A estatística descritiva e inferencial pode ser obtida com ANOVA1f_indep_FisherWhite_sodio.R:
# ANOVA1f_indep_FisherWhite_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)
library(psych)
library(lattice)
library(car)
library(multcomp)
library(gplots)
library(emmeans)
library(readxl)
source("eiras_plotIC.R")
# suppress warnings
options(warn=-1)
TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))
print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
error.bars="conf.int", level=.95,
connect=FALSE,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")],
main="IC95%",
barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")
cat("\nANOVA unifatorial independente de Fisher",
"\ncom ajuste para heterocedasticidade de White:\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
modelo <- lm(Sodium~Instructor, data=TH)
print(res <- car::Anova(modelo, type=2, white.adjust=TRUE))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- res$F[1]
dfn <- res$Df[1]
dfd <- res$Df[2]
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
"de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
c("H0", "Obs",
paste("Fc(",dfn,",",dfd,") = ",round(fc,3),sep=""),
paste("Fobs = ",round(F,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
cat(paste("Teste post hoc\n"))
print(EMM <- emmeans::emmeans(modelo, "Instructor"))
print(grf <- plot(EMM, colors = "black",
main="Estimated Marginal Means",
xlab=VD,
ylab=VI))
# Testes post hoc
cat("\n\nTestes post hoc:\n")
# nomes do fator encurtados
TH$Instructor <- as.character(TH$Instructor)
fatores <- unique(as.character(TH$Instructor))
letra <- "A"
legenda <- c()
cat ("\nLegenda:\n")
for( f in 1:length(fatores))
{
cat("\t",letra," ... ",fatores[f],"\n",sep="")
legenda <- c(legenda,paste(letra," ... ",fatores[f],"\n",sep=""))
TH$Instructor[TH$Instructor==fatores[f]] <- letra
ascii <- strtoi(charToRaw(letra),16L)
letra <- rawToChar(as.raw(ascii+1))
}
TH$Instructor <- as.factor(TH$Instructor)
modelo <- lm(Sodium~Instructor, data=TH)
mc.tukey <- multcomp::glht(modelo, linfct = mcp(Instructor = "Tukey"))
print(mcs.tukey <- summary(mc.tukey, test=adjusted("bonferroni")))
multcomp::cld(mcs.tukey, level=alfa, decreasing=TRUE)
# plot(mc.tukey,las=3)
t.ic <- mcs.tukey$test$qfunction(0.95)
df_plot <- data.frame(
names (mcs.tukey$test$coefficient),
as.vector(mcs.tukey$test$coefficients),
as.vector(mcs.tukey$test$coefficients)-t.ic*as.vector(mcs.tukey$test$sigma),
as.vector(mcs.tukey$test$coefficients)+t.ic*as.vector(mcs.tukey$test$sigma),
as.vector(mcs.tukey$test$pvalues)
)
eiras_plotIC(df_plot,
main="95% family-wise confidence level",
xlab="Difference",
usecolor="n"
)
legend("topleft",legenda,lwd=0,lty=0,cex=0.6,box.lwd=0, border="transparent", bg="transparent")
mc.dunnett <- multcomp::glht(modelo, linfct = mcp(Instructor = "Dunnett"))
mcs.dunnett <- summary(mc.dunnett, test=adjusted("bonferroni"))
print(mcs.dunnett)
# plot(mc.dunnett,las=3)
t.ic <- mcs.dunnett$test$qfunction(0.95)
df_plot <- data.frame(
names (mcs.dunnett$test$coefficient),
as.vector(mcs.dunnett$test$coefficients),
as.vector(mcs.dunnett$test$coefficients)-t.ic*as.vector(mcs.dunnett$test$sigma),
as.vector(mcs.dunnett$test$coefficients)+t.ic*as.vector(mcs.dunnett$test$sigma),
as.vector(mcs.dunnett$test$pvalues)
)
eiras_plotIC(df_plot,
main="95% family-wise confidence level",
xlab="Difference",
usecolor="n"
)
legend("topleft",legenda,lwd=0,lty=0,cex=0.6,box.lwd=0, border="transparent", bg="transparent")
# enable warnings
options(warn=0)Obtendo:
Loading required package: carData
Attaching package: 'car'The following object is masked from 'package:psych':
logitLoading required package: mvtnormLoading required package: survivalLoading required package: TH.dataLoading required package: MASS
Attaching package: 'TH.data'The following object is masked from 'package:MASS':
geyser
Attaching package: 'gplots'The following object is masked from 'package:stats':
lowess
Descriptive statistics by group
group: Brendon Small
vars n mean sd median trimmed mad min max range skew
X1 1 20 1287.5 193.73 1300 1284.38 166.79 950 1700 750 0.12
kurtosis se
X1 -0.46 43.32
--------------------------------------------------------
group: Coach McGuirk
vars n mean sd median trimmed mad min max range skew
X1 1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525 525 0.3
kurtosis se
X1 -0.85 31.84
--------------------------------------------------------
group: Melissa Robins
vars n mean sd median trimmed mad min max range skew
X1 1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400 500 0.08
kurtosis se
X1 -1.1 32.01
ANOVA unifatorial independente de Fisher
com ajuste para heterocedasticidade de White:
VD = Sodium
Fator = Instructor
Analise de significancia estatistica: teste omnibusCoefficient covariances computed by hccm()
Analysis of Deviance Table (Type II tests)
Response: Sodium
Df F Pr(>F)
Instructor 2 5.5748 0.006148 **
Residuals 57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.1636056
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor
Selecao de modelo
- R^2 ajustado = 0.0731098
- omega^2 = 0.07195982
Teste post hoc
Instructor emmean SE df lower.CL upper.CL
Brendon Small 1288 36.1 57 1215 1360
Coach McGuirk 1246 36.1 57 1174 1319
Melissa Robins 1124 36.1 57 1051 1196
Confidence level used: 0.95 Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
Testes post hoc:
Legenda:
A ... Brendon Small
B ... Coach McGuirk
C ... Melissa Robins
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Sodium ~ Instructor, data = TH)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
B - A == 0 -41.25 51.09 -0.807 1.00000
C - A == 0 -163.75 51.09 -3.205 0.00664 **
C - B == 0 -122.50 51.09 -2.398 0.05939 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: lm(formula = Sodium ~ Instructor, data = TH)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
B - A == 0 -41.25 51.09 -0.807 0.84559
C - A == 0 -163.75 51.09 -3.205 0.00443 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)
Analise a saída deste RScript e encontre:
- a estatística descritiva (textual e gráficos)
- a ANOVA unifatorial
- encontre os graus de liberdade \(F\)(2, 57).
- o valor de \(F\) calculado, igual a 5.5748, NA.
- o valor de \(p\) correspondente, igual a 0.006148, NA.
- o tamanho de efeito, \(\eta^2\)=0.1636056.
- os testes post hoc com a legenda:
A ... Brendon Small
B ... Coach McGuirk
C ... Melissa Robins
- e seus resultados em:
- contrastes de Tukey (textual e gráfico) com os respectivos valores \(p\): 1, 0.0066393, 0.0593866
- contrastes de Dunnett (textual e gráfico) com os respectivos valores \(p\): 0.8455883, 0.0044262
- note que \(p=1\) não pode existir, é um arredondamento; por este motivo, no gráfico, fizemos uma correção.
|
Existem muitas famílias de testes para análise post hoc disponíveis. Uma família de testes para análise post-hoc mais conservadora é adequada quando há um grande número de comparações como, por exemplo, o das diferenças honestamente significantes (DHS) de Tukey. “O Teste proposto por Tukey (1953) é também conhecido como teste de Tukey da diferença honestamente significativa (honestly significant difference)(HSD) e teste de Tukey da diferença totalmente significativa (wholly significant difference)(WSD). É um teste exato em que, para a família de todas as \(c=\frac{1}{2}k(k-1)\) comparações duas a duas, a taxa de erro da família dos testes (FWER) é exatamente \(\alpha\) (e o intervalo de confiança é exatamente \(1-\alpha\)). Métodos de comparações múltiplas exatos são raros. O teste de Tukey tem sido mostrado analiticamente ótimo, no sentido que, entre todos os procedimentos que resultam em intervalos de confiança com mesmo tamanho para todas diferenças duas a duas com coeficiente de confiança da família de pelo menos \(1-\alpha\), o teste de Tukey resulta em intervalos menores. Isso quer dizer que, se a família consiste em todas comparações duas a duas e o teste de Tukey pode ser usado, ele resultará em intervalos menores que qualquer outro método de comparação múltipla de uma etapa.” |
Exemplo de relatório para a ANOVA unifatorial independente
Há 60 participantes no estudo balanceado com três condições independentes, sendo que nenhuma delas é de controle. Os participantes foram distribuídos aleatoriamente e balanceadamente nos três grupos. As suposições de normalidade e homocedasticidade da VD foram consideradas válidas. As médias amostrais brutas mostram que menor ingesta de sódio foi observada entre estudantes submetidos ao programa aplicado por Melissa Robins. A ingesta de sódio dos estudantes de Brendon Small e do Coach McGuirk são semelhantes. A análise de variância de um fator fixo entre participantes mostrou que o efeito fixo Sódio é estatisticamente significante, pois o teste omnibus produziu F(2;57) = 5,57 e p = 0,00615. O tamanho do efeito do fator Instrutor é expresso por eta ao quadrado, sendo que seu valor é igual a 0,16. Portanto, 16% da variância da ingesta de sódio é explicada pelo programa adotado pelos respectivos instrutores. Os testes post hoc de Tukey-HSD confirmaram que as diferenças entre os programas adotados por Robins e Small, e entre os adotados por Robins e McGuirk são estatisticamente significantes. Não se observou diferença estatisticamente significante entre os programas adotados por Small e McGuirk.
ANOVA unifatorial independente de Welch
O RScript ANOVA1f_indep_Welch_sodio.R é uma modificação mais robusta do anterior.
As linhas de código que produzem a ANOVA de Fisher-White:
modelo <- lm(Sodium~Instructor, data=TH)
print(res <- car::Anova(modelo, type=2, white.adjust=TRUE))são substituídas por:
print(res <- jmv::anovaOneW(formula = Sodium~Instructor, data=TH,
desc=TRUE, descPlot = FALSE, phMethod ='gamesHowell',
phMeanDif = TRUE, phTest=TRUE, phFlag=TRUE))O objeto res, que recebe o que retorna da função jmv::anovaOneW() é diferente de car::Anova() e, portanto, na sequência do RScript, a captação dos valores foram devidamente ajustados. Esta função da library jmv já executa os testes post hoc.
Executando:
# ANOVA1f_indep_Welch_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)
library(psych)
library(lattice)
library(car)
library(gplots)
library(jmv)
library(readxl)
# suppress warnings
options(warn=-1)
TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))
print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
error.bars="conf.int", level=.95,
connect=FALSE,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")],
main="IC95%",
barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")
cat("\nANOVA unifatorial independente de Welch\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
print(res <- jmv::anovaOneW(formula = Sodium~Instructor, data=TH,
desc=TRUE, descPlot = FALSE, phMethod ='gamesHowell',
phMeanDif = TRUE, phTest=TRUE, phFlag=TRUE))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- as.numeric(res$anova$asDF[2])
dfn <- as.numeric(res$anova$asDF[3])
dfd <- as.numeric(res$anova$asDF[4])
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
"de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
c("H0", "Obs",
paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""),
paste("Fobs = ",round(F,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)Obtém-se:
Attaching package: 'jmv'The following object is masked from 'package:psych':
pcaThe following object is masked from 'package:stats':
anova
Descriptive statistics by group
group: Brendon Small
vars n mean sd median trimmed mad min max range skew
X1 1 20 1287.5 193.73 1300 1284.38 166.79 950 1700 750 0.12
kurtosis se
X1 -0.46 43.32
--------------------------------------------------------
group: Coach McGuirk
vars n mean sd median trimmed mad min max range skew
X1 1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525 525 0.3
kurtosis se
X1 -0.85 31.84
--------------------------------------------------------
group: Melissa Robins
vars n mean sd median trimmed mad min max range skew
X1 1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400 500 0.08
kurtosis se
X1 -1.1 32.01
ANOVA unifatorial independente de Welch
VD = Sodium
Fator = Instructor
Analise de significancia estatistica:
One-Way ANOVA (Welch's)
dep F[welch] df1[welch] df2[welch] p[welch]
"Sodium" Sodium 5.765465 2 37.39637 0.006569226
Group Descriptives
dep group num mean sd se
"SodiumBrendon Small" Sodium Brendon Small 20 1287.50 193.7341 43.32026
"SodiumCoach McGuirk" Sodium Coach McGuirk 20 1246.25 142.4123 31.84435
"SodiumMelissa Robins" Sodium Melissa Robins 20 1123.75 143.1495 32.00920
Games-Howell Post-Hoc Test – Sodium
Brendon Small Coach McGuirk Melissa Robins
Brendon Small Mean difference - 41.25 163.75
t-value - 0.7672 3.0401
df - 34.8931 34.9827
p-value - 0.7253 0.0121
Coach McGuirk Mean difference - - 122.5
t-value - - 2.7131
df - - 37.999
p-value - - 0.0263
Melissa Robins Mean difference - - -
t-value - - -
df - - -
p-value - - -
Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.2356748
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 
Selecao de modelo
- R^2 ajustado = 0.1104092
- omega^2 = 0.1079091 Analise a saída deste RScript e encontre:
- a mesma estatística descritiva do RScript anterior.
- os valores da ANOVA de Welch, com o valor \(F\text{ calculado}\)=(2,37.4)=5.7655 e o valor-\(p\)=0.0065692 (note o número fracionário para os graus de liberdade do denominados, da mesma forma que acontece com o teste \(t\) de Welch).
- os testes post hoc, na forma de uma matriz relacionando os pares de programas (Small, McGuirk e Robins), mostrando:
- a diferença média
- o valor da estatística \(t\) (são feitos testes \(t\) par-a-par)
- os graus de liberdade fracionários (\(t\) de Welch)
- os valores \(p\), já corrigidos para serem interpretados em comparação com \(\alpha\)=0.05.
- o valor estimado para o tamanho de efeito do teste omnibus, dado por \(\eta^2\)=0.2356748.
- note que este procedimento não fornece os IC95% pós-modelo; pretendemos incorporar isto depois
A conclusão, neste exemplo, é similar à obtida com o teste de ANOVA de Fisher-White: o programa adotado por Melissa Robins obteve ingestas de sódio significantemente menores que a dos outros dois instrutores.
ANOVA unifatorial por bootstrapping
O RScript ANOVA1f_indep_Bootstrap_sodio.R é a versão com bootstrapping do pacote lmboot.
O RScript que implementa ANOVA desta maneira é:
# ANOVA1f_indep_Bootstrap_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)
library(psych)
library(lattice)
library(car)
library(gplots)
library(lmboot)
library(readxl)
# suppress warnings
options(warn=-1)
TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))
print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
error.bars="conf.int", level=.95,
connect=FALSE,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")],
main="IC95%",
barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")
cat("\nANOVA unifatorial por bootstrapping\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
# remova seed para ter resultados variados
bootsamples <- 1e6
modelo_boot <- lmboot::ANOVA.boot(Sodium~Instructor,
B = bootsamples, type = "residual",
wild.dist = "normal",
# seed = 123,
data = TH,
keep.boot.resp = FALSE)
cat(paste("F(",dfn,",",round(dfd,3),") = ",round(F,5),", p = ",round(modelo_boot$`p-values`,5),"\n", sep=""))
cat(paste("(",bootsamples," bootstrap samples)\n", sep=""))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- qf(1-modelo_boot$`p-values`, modelo_boot$df[1], modelo_boot$df[2])
dfn <- as.numeric(modelo_boot$df[1])
dfd <- as.numeric(modelo_boot$df[2])
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
"de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- modelo_boot$`p-values`
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
c("H0", "Obs",
paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""),
paste("Fobs = ",round(F,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)Obtendo-se:
Descriptive statistics by group
group: Brendon Small
vars n mean sd median trimmed mad min max range skew
X1 1 20 1287.5 193.73 1300 1284.38 166.79 950 1700 750 0.12
kurtosis se
X1 -0.46 43.32
--------------------------------------------------------
group: Coach McGuirk
vars n mean sd median trimmed mad min max range skew
X1 1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525 525 0.3
kurtosis se
X1 -0.85 31.84
--------------------------------------------------------
group: Melissa Robins
vars n mean sd median trimmed mad min max range skew
X1 1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400 500 0.08
kurtosis se
X1 -1.1 32.01
ANOVA unifatorial por bootstrapping
VD = Sodium
Fator = Instructor
Analise de significancia estatistica: teste omnibus
F(2,37.396) = 5.76546, p = 0.00635
(1e+06 bootstrap samples)
Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.1626598
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 
Selecao de modelo
- R^2 ajustado = 0.07253934
- omega^2 = 0.07139764 Verifique na saída:
- compare com as duas soluções anteriores.
- note que este procedimento não fornece os IC95% pós-modelo; pretendemos incorporar isto depois
ANOVA independente, sem os dados brutos
É muito comum, em publicações, que somente tenhamos acesso às medidas-resumo (número de participantes, média, desvio-padrão e correlação). Nestes casos, os RScripts acima não são utilizáveis.
Para fazer os testes quando os dados brutos não estão disponíveis, criamos os seguintes scripts:
ANOVA unifatorial relacionada ou para medidas repetidas
Utilizada quando os participantes são avaliados sob todas as condições experimentais, i.e., um delineamento intraparticipantes.
Neste delineamento a variação entre os grupos não é devida às diferenças individuais (cada participante em cada grupo é o mesmo), e a fórmula para o cálculo da estatística \(F\) leva este fato em consideração.
O delineamento com medidas repetidas leva a um teste com maior poder.
No entanto, como o mesmo participante é controle de si mesmo e submetido a todas as condições experimentais, sugere-se que esta ordem de submissão, se possível, seja aleatorizada.
- suposições
– as diferenças dos valores das VDs são independentes entre as unidades observacionais. – as diferenças dos valores das VDs têm distribuição normal multivariada. – esfericidade: homocedasticidade das variâncias das diferenças das VDs (uma explicação sem muita matemática está em https://en.wikipedia.org/wiki/Mauchly%27s_sphericity_test)
“O teste t é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que [para um delineamento] completamente intraparticipantes, [deve haver] pelo menos 22 participantes ao todo.” (Dancey & Reidy, 2019, p. 472)
- situação
| Usaremos os mesmos dados do exemplo anterior, mas imaginando que as 20 medidas de ingesta de sódio sejam do mesmo participante, submetido aos três diferentes programas educacionais. Então, a normalidade multivariada das três diferenças não pode ser assumida, pois há 20 unidades observacionais no estudo. |
O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.
O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.
Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.
- planejamento
Brendon Small, Coach McGuirk e Melissa Robins fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Estudantes atenderam os diferentes programas de educação nutricional suscessivamente, e os instrutores querem ver se a ingestão média de sódio é a mesma quando cada um dos três programas foi seguido.
- hipóteses nula e alternativa
O delineamento do estudo é diferente, mas as hipóteses são as mesmas:
\[H_0: \mu_{\text{Small}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\] \[H_1: \exists |\mu_i \ne \mu_j; i \ne j; i,j=1,2,3\]
Operacionalização do teste
ANOVA unifatorial relacionada ou para medidas repetidas, balanceadas
A planilha a ser utilizada é Nutricao3par.xlsx. Verifique a coluna Student, que foi alterada para indicar o participante nos três programas (compare com Nutricao3.xlsx, usada para ANOVA independente).
A libraries necessárias são lmerTest e ez e o trecho que aplica o método estatístico é:
res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)Executando:
# ANOVA1f_dep_balanc_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)
library(psych)
library(lattice)
library(car)
library(lmerTest)
library(ez)
library(readxl)
# suppress warnings
options(warn=-1)
TH <- read_excel("Nutricao3par.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))
print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
jitter.x=TRUE, col="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")
cat(paste("\nTeste omnibus por GLMM com efeito aleatorio\n"))
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)
cat("\n"); print(res1 <- stats::anova(res)); cat("\n")
print(lmerTest::rand(res))
cat(paste("\nTeste post hoc\n\n"))
print(out <- lmerTest::difflsmeans(res))
dfn <- as.numeric(res1[3])
dfd <- as.numeric(res1[4])
F <- as.numeric(res1[5])
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("\neta^2 =", eta2, "\nGrau", geta2,
"de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
c("H0", "Obs",
paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""),
paste("Fobs = ",round(F,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
cat("Teste omnibus por GLM univariado\n")
# ezANOVA: apenas para medidas repetidas balanceadas
print(res <- ez::ezANOVA(data=TH, dv=Sodium, wid=Student, within=Instructor,
detailed=TRUE, type=3))
# enable warnings
options(warn=0)Obtém-se:
Loading required package: lme4Loading required package: MatrixRegistered S3 methods overwritten by 'lme4':
method from
cooks.distance.influence.merMod car
influence.merMod car
dfbeta.influence.merMod car
dfbetas.influence.merMod car
Attaching package: 'lmerTest'The following object is masked from 'package:lme4':
lmerThe following object is masked from 'package:stats':
step
Descriptive statistics by group
group: Brendon Small
vars n mean sd median trimmed mad min max range skew
X1 1 20 1287.5 193.73 1300 1284.38 166.79 950 1700 750 0.12
kurtosis se
X1 -0.46 43.32
--------------------------------------------------------
group: Coach McGuirk
vars n mean sd median trimmed mad min max range skew
X1 1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525 525 0.3
kurtosis se
X1 -0.85 31.84
--------------------------------------------------------
group: Melissa Robins
vars n mean sd median trimmed mad min max range skew
X1 1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400 500 0.08
kurtosis se
X1 -1.1 32.01
Teste omnibus por GLMM com efeito aleatorio
VD = Sodium
Fator = Instructor
Analise de significancia estatistica: teste omnibus
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 290146 145073 2 38 30.581 1.217e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA-like table for random-effects: Single term deletions
Model:
Sodium ~ Instructor + (1 | Student)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 5 -352.02 714.05
(1 | Student) 4 -375.21 758.42 46.376 1 9.763e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Teste post hoc
Least Squares Means table:
Estimate Std. Error df
InstructorBrendon Small - InstructorCoach McGuirk 41.2500 21.7807 38
InstructorBrendon Small - InstructorMelissa Robins 163.7500 21.7807 38
InstructorCoach McGuirk - InstructorMelissa Robins 122.5000 21.7807 38
t value lower
InstructorBrendon Small - InstructorCoach McGuirk 1.8939 -2.8426
InstructorBrendon Small - InstructorMelissa Robins 7.5181 119.6574
InstructorCoach McGuirk - InstructorMelissa Robins 5.6243 78.4074
upper Pr(>|t|)
InstructorBrendon Small - InstructorCoach McGuirk 85.3426 0.06587 .
InstructorBrendon Small - InstructorMelissa Robins 207.8426 4.949e-09 ***
InstructorCoach McGuirk - InstructorMelissa Robins 166.5926 1.866e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
Degrees of freedom method: Satterthwaite
Analise de significancia pratica: tamanho de efeito
eta^2 = 0.6167848
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 
Selecao de modelo
- R^2 ajustado = 0.4313252
- omega^2 = 0.4251262
Teste omnibus por GLM univariado
$ANOVA
Effect DFn DFd SSn SSd F p p<.05
1 (Intercept) 1 19 89182041.7 1307541.7 1295.91189 6.002751e-19 *
2 Instructor 2 38 290145.8 180270.8 30.58049 1.217342e-08 *
ges
1 0.9835909
2 0.1631905
$`Mauchly's Test for Sphericity`
Effect W p p<.05
2 Instructor 0.3889714 0.0002038255 *
$`Sphericity Corrections`
Effect GGe p[GG] p[GG]<.05 HFe p[HF]
2 Instructor 0.6207214 3.881854e-06 * 0.6427565 2.772133e-06
p[HF]<.05
2 *Analise a saída deste RScript e encontre:
- a estatística descritiva (sempre a mesma).
- observe o valor de \(F\text{ calculado}\) e valor-\(p\) correspondentes (neste exemplo fictício, os dados são os mesmos da ANOVA independente).
- localize e interprete os testes post hoc; compare com os valores-\(p\) da ANOVA independente.
- aparece um novo trecho, “Teste omnibus por GLM univariado”.
ANOVA unifatorial relacionada ou para medidas repetidas, desbalanceadas
Vamos usar os mesmos dados, mas eliminar uma única medida (removemos o estudante a de Melissa Robins, planilha Nutricao3parD.xlsx).
Para lidar com esta situação, o RScript execute:
# ANOVA1f_dep_desbalanc_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)
library(psych)
library(lattice)
library(car)
library(lmerTest)
library(readxl)
# suppress warnings
options(warn=-1)
TH <- read_excel("Nutricao3parD.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))
print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
ylab=names(TH)[which(names(TH)=="Sodium")],
xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
jitter.x=TRUE, col="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")
cat(paste("\nTeste omnibus por GLMM com efeito aleatorio\n"))
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)
cat("\n"); print(res1 <- stats::anova(res)); cat("\n")
print(lmerTest::rand(res))
cat(paste("\nTeste post hoc\n\n"))
print(out <- lmerTest::difflsmeans(res))
dfn <- as.numeric(res1[3])
dfd <- as.numeric(res1[4])
F <- as.numeric(res1[5])
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("\neta^2 =", eta2, "\nGrau", geta2,
"de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
c("H0", "Obs",
paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""),
paste("Fobs = ",round(F,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)É o mesmo RScript usado para a ANOVA balanceada, exceto pelas linhas finais (“Teste omnibus por GLM univariado”). Obtém-se:
Descriptive statistics by group
group: Brendon Small
vars n mean sd median trimmed mad min max range skew
X1 1 20 1287.5 193.73 1300 1284.38 166.79 950 1700 750 0.12
kurtosis se
X1 -0.46 43.32
--------------------------------------------------------
group: Coach McGuirk
vars n mean sd median trimmed mad min max range skew
X1 1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525 525 0.3
kurtosis se
X1 -0.85 31.84
--------------------------------------------------------
group: Melissa Robins
vars n mean sd median trimmed mad min max range skew
X1 1 19 1135.53 136.76 1125 1132.35 148.26 925 1400 475 0.08
kurtosis se
X1 -1.04 31.37
Teste omnibus por GLMM com efeito aleatorio
VD = Sodium
Fator = Instructor
Analise de significancia estatistica: teste omnibus
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 254746 127373 2 37.086 27.369 4.998e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA-like table for random-effects: Single term deletions
Model:
Sodium ~ Instructor + (1 | Student)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 5 -345.58 701.16
(1 | Student) 4 -368.17 744.34 45.177 1 1.8e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Teste post hoc
Least Squares Means table:
Estimate Std. Error
InstructorBrendon Small - InstructorCoach McGuirk 41.2500 21.5730
InstructorBrendon Small - InstructorMelissa Robins 157.5073 21.9807
InstructorCoach McGuirk - InstructorMelissa Robins 116.2573 21.9807
df t value lower
InstructorBrendon Small - InstructorCoach McGuirk 37.0 1.9121 -2.4601
InstructorBrendon Small - InstructorMelissa Robins 37.1 7.1657 112.9749
InstructorCoach McGuirk - InstructorMelissa Robins 37.1 5.2891 71.7249
upper Pr(>|t|)
InstructorBrendon Small - InstructorCoach McGuirk 84.9601 0.06362 .
InstructorBrendon Small - InstructorMelissa Robins 202.0397 1.677e-08 ***
InstructorCoach McGuirk - InstructorMelissa Robins 160.7897 5.709e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
Degrees of freedom method: Satterthwaite
Analise de significancia pratica: tamanho de efeito
eta^2 = 0.5961188
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 
Selecao de modelo
- R^2 ajustado = 0.4091071
- omega^2 = 0.4028569 Exemplo de relatório para a ANOVA unifatorial relacionada
Vinte participantes foram selecionados para um estudo com delineamento intrapartipantes com três condições experimentais. A ordem de aplicação das três condições experimentais em cada participante foi aleatorizada. As médias amostrais brutas mostram que menor ingesta de sódio foi observada entre estudantes submetidos ao programa aplicado por Melissa Robins. A ingesta de sódio dos estudantes de Brendon Small e do Coach McGuirk são semelhantes. A ANOVA unifatorial para medidas repetidas mostrou que qualquer diferença entre as condições experimentais é improvável de ter ocorrido apenas por erro amostral, considerando a hipótese nula verdadeira, pois a estatística de teste multivariado observada é F(2;38) = 30,581 e o valor-p associado é igual a \(1.22 \cdot 10^{-8}\). O tamanho do efeito do fator intraparticipantes álcool é estimado pelo eta ao quadrado cujo valor indica que 61,7% da variância da ingesta de sódio é explicada pelo efeito do fator fixo Instrutor. A análise post-hoc confirmou que as diferenças entre entre os programas adotados por Robins e Small, e entre os adotados por Robins e McGuirk são estatisticamente significantes. Não se observou diferença estatisticamente significante entre os programas adotados por Small e McGuirk.